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Abstract 

We present analytical and numerical treatments for evaluating the 
time-of-flight momentum distribution for the stationary states of a two- 
boson system trapped in a quartic double- well potential, paying particular 
attention to the Tonks and noninteracting regimes. We find that the time- 
of-flight distributions can serve as a valuable tool in profiling the states of 
this system, and provide the tof plots for low-energy excitations using a 
nonadaptive finite element method that is more efficient than traditional 
finite difference methods. 

1 Introduction 

The study of ultracold boson systems in Optical Dipole Traps (ODT) such as 
Bose-Einstein condensates, have received considerable attention in academia 
over the last decade of the 20*'' century [1] [2] [3] [4]. More recently, experi- 
ments involving number squeezed states of trapped alkali atoms have yielded 
promise [5]. New techniques, such as quantum tw^eezing [6] and quantum many- 
body culling [5] , are being developed that can create mesoscopic two-boson sys- 
tems out of ultracold atoms in optical traps. Theoretical studies demonstrate 
the possibility of number state generation by atomic culling as well, where a 
BEC is number-squeezed by "culling" atoms from a trapped condensate down 
to a sub-poissonian regime, making the number uncertainty small enough to be 
ignored [5]. Such a two-boson system can be subjected to a micrometer-scale 
double well by various means, ranging from small volume optical traps [5], to 
atom chips [7] [8] . An optical lattice of such double- wells can also be generated 
by two counter-propagating lasers of linearly polarized light with a known an- 
gle between their planes of polarization, and a transverse magnetic field to mix 
the two potentials [9]. If the on-site lattice depth is sufficiently deep then the 
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tunneling between the sites can be neglected. Furthermore, if they are loaded 
homogeneously from a cold-atom system confined in an optical dipole trap by 
atom culling [5], each double well system can be treated in isolation exactly 
as depicted in [10]. More recently (2009), number squeezing and subpoissonian 
distribution of atoms in each site in an optical lattice have been reported by Itah 
et al [11]. The weak nature of the interactions of such atoms means that such 
cold atom systems are useful tools for quantum information processing [12], as 
well as in studying quantum entanglement [13]. More recently, the study of the 
dynamics of quantum control in such systems have been performed as well [10]. 
The momentum distributions of two-bosons systems have been evaluated nu- 
merically by Murphy and McCann using a finite difference method [14]. We 
present the numerical evaluation of the time of flight momentum distributions 
of the lower energy stationary states of this system using a nonadaptive finite 
element method that is conputationally faster and more accurate. 

In the following sections, we evaluate the time of flight (tof) signatures of 
these wavefunction. Section [5] details how the double well system was diago- 
nalized and the eigenfunctions obtained. Section [3] discusses the nature of the 
time-of-flight signatures of the different states, and numerical results are shown 
in section [D Concluding remarks are made in the final section. 

2 The Eigensystem: Strongly Interacting and 
Single Particle Regimes 

Our system consists of two alkali metal bosons confined to a double-well optical 
potential. The effective interaction between the bosons, in three dimensions, is 
obtained in the long wavelength approximation to be 



where h is Planck's constant, as is the s-wave scattering length and = 
{xi,yi,Zi) is the displacement of the ith particle [15] [16]. The system can 
be confined in two spatial (radial) directions so that the essential dynamics oc- 
curs in the x - direction by the use of anisotropic magnetic traps with high 
aspect ratio [17] [18]. In that case, the other 2 dimensions can be integrated 
out [10] [17], yielding an effective 1-dimensional interaction 



We will consider the case of two identical bosons confined to a quartic double 
well potential. We get the total Hamiltonian for the system to be 



H =pI+pI + Voi-2xl + xj) + Voi-2xl + x^) + Uo6{xi - X2). (3) 



where pi is the momentum of the ith particle {i— 1,2), Xi is the position of 
the ith particle along the x-axis, and Vq determines the depth of the double well 
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potential. We have used dimensionless expressions for all the degrees of freedom, 
as well as the system parameters, by introducing a characteristic length scale 
Lu- Thus, the actual Haniiltonian H' relates to the dimensionless Hamiltonian 
H as H = where £'„ — 2mL'^ ■ Similarly, Ua — "^"-^""^ and the time scales 



Fig. [T] shows a plot of the quartic double well 



as t — where — 
V{x) = Voi-'2x^ + x^) for weU depth Vq = 7.2912229. 

The numerical diagonalization of the Hamiltonian in Eq. [3] is facilitated by 
a nonadaptive finite element method using the analytically obtained matrix 
elements of the Hamiltonian in a finite wave train basis of size L — 3.5 (in units 
of Lu), 
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[{xi\ni){x2\n2) + {xi\n2){x2\ni)]. (4) 



Here, 



Vl 



sm 



mr X 



1) 



(5) 



within the range —L < x < L and vanishes outside. 

We will investigate the tof distributions in two regimes of the {Vq, Uq) pa- 
rameter space of the double well system. Here, Vq is the well depth, and Uq 
the amplitude of the point contact pseudopotential in 1-dimension. The first 
regime, henceforth referred to as the 'strongly interacting regime' will consist 
of a very strongly repulsive system and a moderate well depth. We define the 
'strongly interacting factor' for this system, 7, as 



Uo 



(6) 



Here, E, the energy of the state, is a measure of the ability of the bosons to 
tunnel across from one well to another. When 7 — > 00, we reach the strongly 
interacting regime where the interaction completely dominates the system [19]. 
Figure [5] shows the evolution of the ground state of the system as 7 is increased. 
The order parameter being plotted as a function of 7 for the ground state is 
Pi/ 51, where 

p, = 8l [ dx\{x,x\Ei)\\ (7) 



Here, i is an index distinguishing different regimes of interest in the 7-space. 
Also, Pi is the total probability that the two particles will be together within 
a rectangular strip along the line xi — X2 and arbitrarily small width 51. As 
expected, it vanishes for large values of 7. 

In this strongly interacting regime, the two particles have no probability of 
occupying the same position simultaneously. Thus, they act in a way that is 
similar to a Tonks gas [19]. The transition to this regime is not consistent, 
however. We note four distinct ranges of 7 for which the decay rates of pi/5l 
are different. In the first three ranges, pi seems to be decaying exponentially 
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Figure 1: Plot of the double- well potential experienced by each boson for the 
single particle regime. The energy levels, Ei — —6.42262, E2 = —5.68883 
and E4 = 0.640055 of the interacting two-boson system (interaction strength 
Uo = -1.0) are also sketched. Here, Vq = 7.2912229. 
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Figure 2: Semi-logarithmic plot of the one dimensional probability density Pi/5l 
of two particles being within a rectangular strip of arbitrarily small width 51 
along the line x\ = X2- Pi/Sl is plotted as a function of the strongly interacting 
parameter 7 for a constant Vq = 4.0. The decay rate of the probability pi 
changes sharply at 4 regions, labeled by the index i. The data points (indicated 
by circles) have been fitted to exponential decay rates at each region (indicated 
by lines). The legend provides the numerically fitted values of the decay rates 
7?. Note the discontinuous spike at 7 7. 
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ie (pi/Sl) = ip^/Sl) e~^^/''^'' for i = 1,2,3. The data points have been fitted 
to exponents by the use of numerical nonhnear least-squares algorithms. The 
decay rate, characterized by 7°, decreases sharply at 7 ~ 1,2 and 6. Near 
7 ~ 7, there is a sharp increase in pi after which it continues to decrease. If we 
neglect the probability if it falls below 1/e of the maximum, then the 'strongly 
interacting regime' is achieved beyond 7 ~ 0.4. In our case, we have chosen a 
7 of 5.20142 for our strongly interacting regime, placing the system in region 3 
of FigH The value of (Vo, Uq) chosen is (4.0,40.0). 

The second regime, henceforth referred to as the 'single particle regime', will 
consist of a weakly attractive system and the well-depth as seen in [10]. Thus, 
the parameter values chosen are (7.2912229, —1.0). The probability distributions 
of the ground state as well as the excited states |i?2) and |£'4), given by 
Eqn [3] are shown in Figs [Sja through [3]c for the strongly interacting regime. 
Note that, as expected, there is virtually no probability that xi = X2- The 
probability distributions of the first seven quantum energy states of the system 
in the single particle regime are shown in Figs. [Ha through [Hg. Note the plots 
of the ground state, |i?i), third excited state, I-E4), and sixth excited state \E^). 



The normalized first order correlation function of a single double well is a mea- 
surement of the atomic density n{x). Such correlations can be measured by the 
time-of-flight (TOF) technique in which the trapped atoms are released suffi- 
ciently quickly that the diabatic approximation in quantum mechanics can be 
applied. The atoms then expand ballistically until they reach a detection plate. 
If the plate is far enough from the double well system that the far-field approxi- 
mation can be used, then the Green's Function for the system can be simplified 
and the time translation reduced to a simple Fourier Transform. The 'detector 
plane' coordinates are denoted by unprimed variables [xi,X2,t] and the double 
-well coordinates by primed variables [x'l , X2 , t] for a 2-particle system after all 
the external fields and traps have been diabatically switched off. The interac- 
tions between the atoms while in flight can be rendered negligible by tuning a 
homogeneous magnetic field close to the Feshbac resonance that adds an attrac- 
tive amplitude to the normally repulsive point contact pseudopotential [20] [16]. 
The system then evolves ballistically in free space. 

The Green's Function or Propagator (^(x, i; x', t') is defined by 



where x = [xi, X2\, and ^'(x, t) is the wavefunction, with similar expressions for 
the primed coordinates. For free space, the relevant 1-dimensional Schrodinger 



3 Time of Flight Images 
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Figure 3: Plots of energy eigenfunctions for the two interacting bosons in a 
double well potential in the strongly interacting regime. Figures (a) through 
(c) are contour plots of the probability density \{xi, X2\Ei)\'^ , \{xi, X2\E2)\'^ and 
|(xi,X2|-E4)P respectively. The probabilities are plotted as functions of xi and 
X2- AH units for all figures are dimcnsionless 
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Figure 4: Plots of energy eigenfunctions for the two interacting bosons in a 
double well potential in the single particle regime. All units are dimensionless. 
Figures (a) through (f) are contour plots of the probability density | (xi , a;2 |£^i) ^ 
through \{xi,X2\Eq)\'^ respectively. Figure (g) is a contour plot of the proba- 
bility density \{xi,X2\Ej)\'^.Th.e peaks in the probability are numbered. The 
probabilities are plotted as functions of xi and X2. All units for all figures are 
dimensionless 
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equation for 2-particles is 



H 



H 



dt 



52 



,2 n 



Thus, the Green's function [21] will be the solution to 

d' 



H 



dt 



G(x, t; x', t') = S{x - x')S{t - t'). 



The solution to Eqn[Tn]in free space is 



G(x,i;x',t') 



L 



■ exp 



Here, the ballistic de-Broglie equation. 



= 47rr, 



(9) 



(10) 



(11) 



(12) 



provides the relation between the detector-system separation L and the time- 
of-flight T = {t-t'). 

Now, consider such a two particle system localized at a site j. The wavefunc- 
tion is localized about x^ = [x'j, x'j] and can be written in the form 5'(x' — x^ ). 
We now use Eqns [TT] and [T^ on Eqn [HI and apply the shift theorem for Fourier 
transforms [22] [23] to get 



2t 



ixr 

2 



i^[*(x')]„. . , 



(13) 



where the primed coordinates refer to the double well system, the unprimed 
coordinates refer to the detector, and F [^'(x')]^ is the Fourier transform 



x')e™'^'. 



(14) 



In the equation above, u = [ui, U2] is the momentum space vector. For a large 
collection of such systems, each in the desired pure state, the measured TOF is 
simply the probability obtained from Eqn [13] times the number of such double 
wells N (which we shall subsequently drop off as an appropriately adjusted 
overall normalization) . 



ni^)=N—\F[^i^')]^^^ 



(15) 



In the next section, Eqn [T5| will be evaluated numerically for two bosons in 
a double well, and the density functional 



n(x)^Jdx"ni^,x") (16) 
will be calculated in order to provide the 1-dimensional tof distribution. 
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4 Time-of-Flight: Numerical Plots 

This section will detail the procedure for obtaining numerical plots of the tof 
distributions of the eigenstates of the double well system. The two boson prob- 
lem in a double well is diagonalized as detailed in section [2] and [10]. Thus, the 
eigenfunctions are obtained as linear superpositions of the eigenfunctions of two 
bosons in a box of appropriately chosen length L, ie 



{xi,X2\ni,n2)'"''^ 



1 



^/2(1 + 5ni,n2) 



(a;i|ni)(a;2|n2) + {xi\n2) {x2\ni)], 



(x\n) = —= sm 



2 ' 



(17) 



if \xi\ < L, and otherwise. Thus, the final solution to an eigenfunction \En) of 
the double well will be a linear superposition of the 'finite wave train' functions 
defined above, ie 



[N,N] 

{xi,X2\En)^ ^ C^E,^'"'^\xi,X2\ni,n2)''''\ 
[ni,n2] = [l,l] 



(18) 



where the (7^1 '"^1 are obtained numerically using the nonadaptive finite element 
method. This result can then be substituted into Eqn [15] to get 



n{x,,X2)^N^\ Y: 4", 



ni,n2] 



F 



{x[,x'2\ni,n2) 



[ni,n2] = [l,l] 

Using the linearity of Fourier Transforms and Eqn 1171 we get 

1 



T \x 1 .xo] 



(19) 



F 



{x[,X2\ni,n2Y'^ 



iF[{xi\ni)UF[{x2\n2)U + F[{xi\n2)UF[{x2\ni)]u2) ■ (20) 



The Fourier transform of the finite wave train {{x\n) in EgnllZp can be calculated 
using Gaussian integrations [24] to yield 



F[{x'\n)l 



sin (uL + ^) ,„s sin (uL — ^) ,„s 



(21) 

where u is the momentum space vector. Thus, by plugging Eqn [21] into Eqnl20|, 
and that into Eqn[T9l the tof distribution n{xi,X2) can be obtained, the final 
density distribution is the density functional average of this result viz. 



n{x) — I dx'n{x,x'). 



(22) 
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Thus, a numerical expression for EqnfTSlwas obtained for two degrees of freedom 
xi and X2, and the density functional n{x) determined by integrating out one 
of the coordinates by adaptive Gauss-Kronrod quadrature. 

Numerical results for the tof distributions of the eigenstates of the double 
well for the strongly interacting and single particle regimes are shown in Figs [5] 
and ini respectively. The distributions are shown for tof r — 10'"' units of T„. 
All the dynamics is essentially independent of the characteristic length scale L„ 
(the actual position of the well minima). For practical reasons, we choose an 
Lu of 50 nm [10]. Consequently, with a ®^Rb atomic mass of 85.4678 gmol~^, 
we get a T„ of about 6.7 /zs, which makes r to be 0.67 seconds. Using EgnfT^ 
we get a detector distance of about 2.2 cm. 

Figures [5^ a) through (c) show the tof distributions of the states \Ei), \E2), 
and [£'4) respectively for the strongly interacting gas detailed in section 1. Fig- 
ures |6l[a) through (c) show the tof distributions of the states \Ei), I-E4), and 
lEf) respectively in the single particle regime detailed in section 1. The distru- 
butions match the tof distributions obtained by Murphy and McCann [14], but 
were obtained using a finite element method that is computationally far less 
resource-intensive . 

Time of flight fluorescence methods for profiling the wavefunction, such as 
measuring the momentum distribution by interrupting the particle flow with 
counter-propagating laser beams and then measuring fluorescence as a function 
of time (time of flight absorption) [25] [26] , will have high signal to noise ratio 
(compared to absorption) [5]. Single shot fluorescence images should duphcate 
the profile shown in Figs [6ja)-(c) for a double well system produced by optical 
lattices. For a single magnetically confined double well, repeated measurements 
of position by the means of atom detectors, or by performing scanning tunneling 
microscopy on an appropriate substrate where the atoms are allowed to deposit 
after their tof expansion, should reproduce the required results. 
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5 Conclusions 



We have considered the system of two interacting bosons in a quartic double 
well potential of the type Vq{—2x^+x^). The cigcnspcctrnm has been shown for 
both strongly interacting and nearly single-particle regimes, and the momcntiim 
space distributions shown for both regimes using simple numerical methods. The 
double well arrangement is experimentally obtainable, and the tof distributions 
of the lower energy states are useful information needed for doing quantum 
control and quantum information science problems in such systems. 
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Figure 5: Figures (a) through (c) are plots of the one-dimensional tinic-of-flight 
distributions for the double-well eigenstates \Ei), \E2) and |i?4) respectively in 
the strongly interacting regime. The number density n{x) in the ordinate is for 
10^ double wells after a time of flight r = 10^ (in units of T„). The abscissa is 
shown in dimensionless units of L„. 
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Figure 6: Figures (a) through (c) are plots of the timc-of-flight distributions for 
the double- weU eigenstates |i?4) and \E^) respectively in the single particle 
regime. The number density n{x) in the ordinate is for 10^ double wells after a 
time of flight r = 10^ (in units of r„). The abscissa is shown in dimensionless 
units of L„. 
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